function alphat= filterSmootherCK(data,TT,n,m,k,q,alpha_0,P_0,state_space,varargin)

% State Space system as in Harvey (1989)
% measurement equation: y(t) = d + Z alpha(t) + epsilon(t)
% transition equation: alpha(t+1) = c + T alpha(t) + R eta(t)

alpha_s=zeros(m,TT);

alpha_u=zeros(m,TT+1);
P_u=zeros(m,m,TT+1);
alpha_f=zeros(m,TT);
P_f=zeros(m,m,TT);

%P_transform = @(Pcov) 0.5*(Pcov + Pcov');

alpha_u(:,1)=alpha_0;
P_u(:,:,1)=P_0;



for t=1:TT
   [H,T_star,T,Q_star,Q,Z,d]=state_space(n,k,t,varargin{:});
   
    alpha_f(:,t)=T*alpha_u(:,t);%+c(:,t);
    P_f(:,:,t)=T*P_u(:,:,t)*T'+Q;
    Et=data(:,t)-Z*alpha_f(:,t)-d;
    P_y=Z*P_f(:,:,t)*Z'+H;
    Kt=P_f(:,:,t)*Z'*eye(size(P_y))/(P_y);
    alpha_u(:,t+1)=alpha_f(:,t)+Kt*Et;
    P_u(:,:,t+1)=P_f(:,:,t) - Kt*Z*P_f(:,:,t);
end

alpha_s(:,TT)=mvnrnd1(alpha_u(:,TT+1),P_u(:,:,TT+1),1);

for t=(TT-1):-1:1
    
    Jt=P_u(:,:,t+1)*T_star' * eye(size(Q_star))/(T_star*P_u(:,:,t+1)*T_star'+ Q_star);
    alpha_st=alpha_u(:,t+1)+Jt*(alpha_s(1:k,t+1)-T_star*alpha_u(:,t+1));%-c(1:k,t)
    P_st=P_u(:,:,t+1)-Jt*T_star*P_u(:,:,t+1);
    alpha_s(:,t)=mvnrnd1(alpha_st,P_st,1);
end

ee=[eye(k) zeros(k,k*(q-1))];

alphat=(ee*alpha_s(:,1:end));


